function plotrocks(filenumber)
S=readdata(0,2);
depthrange=[1e3:1e3:4e3 10e3:10e3:100e3]+40e3; 

temphist=zeros(S,length(filenumber),length(depthrange));
refsurface=readdata(0,25); refsurface=refsurface(:,4);
surfacehist=zeros(S,length(filenumber)); avtemphist=surfacehist; timehist=zeros(1,length(filenumber));
timeindex=1;

for count=filenumber
    N=readdata(count,1);
    S=readdata(count,2);
    T=readdata(count,3);
    %     multigridlevels=readdata(count,4);
    %     bound=readdata(count,5);
    time=readdata(count,6);
    %     dtnom=readdata(count,7);
    dx=readdata(count,8);
    dz=readdata(count,9);
    mx=readdata(count,10);
    mz=readdata(count,11);
    %     mcs=readdata(count,12);
    %     mct=readdata(count,13);
    %    mtype=readdata(count,14);
    mTemp=readdata(count,15);
    %mstressxx=readdata(count,16);
    %mstressxz=readdata(count,17);
    %mstrainxx=readdata(count,18);
    %mstrainxz=readdata(count,19);
    %mfinitestrain=readdata(count,20);
    %     vx=readdata(count,21);
    %     vz=readdata(count,22);
    %     P=readdata(count,23);
    Nsurfaces=readdata(count,24);
    surfaces=readdata(count,25);

    count
    x=zeros(1,S); z=zeros(1,T);
    for s=1:S-1 x(s+1)=x(s)+dx(s); end
    for t=1:T-1 z(t+1)=z(t)+dz(t); end
    
    [gridtemp,xgrid,zgrid]=bilingrid(mTemp,mx,mz,dx,dz,1,1);
    [Zgrid,Xgrid] = meshgrid(zgrid,xgrid);
    for d=1:length(depthrange)
        localtemp=interp2(Xgrid',Zgrid',gridtemp',x,depthrange(d));
        temphist(:,timeindex,d)=localtemp;
    end
    integrationdepths=linspace(40e3+0e3,max(z),2*T);
    for s=1:S
        localtemp=interp2(Xgrid',Zgrid',gridtemp',x(s),integrationdepths);
        avtemphist(s,timeindex)=trapz(integrationdepths,localtemp);
    end
    
    sealevel=0.5*(surfaces(1,Nsurfaces)+surfaces(S,Nsurfaces));
    %surfacehist(:,timeindex)=surfaces(:,4);
    %surfacehist(:,timeindex)=surfaces(:,4)-refsurface;
    tempsurface=surfaces(:,4)-0*refsurface; surfacehist(:,timeindex)=tempsurface-mean(tempsurface);  %correction for long term drift
    timehist(timeindex)=time/(60*60*24*365*1e6);
    timeindex=timeindex+1;
end
save surface_time timehist surfacehist temphist avtemphist depthrange x;

load surface_time;
N=length(timehist); S=length(x);
timehistint=linspace(timehist(1),timehist(end),N); dt=timehistint(2)-timehistint(1);

lowcut=0;
for s=1:S
    %surfacehist=interp1(timehist,surfacehist,timehistint); surfacehist(1)=surfacehist(2);
    tempsurface=interp1(timehist,surfacehist(s,:),timehistint); tempsurface(1)=tempsurface(2);
    tempsurface=bandpass(tempsurface,lowcut,1/10,dt);
    surfacehist(s,:)=tempsurface;
    
    for d=1:length(depthrange)
        temptemp=interp1(timehist,temphist(s,:,d),timehistint);
        %temptemp=temptemp-mean(temptemp);
        temptemp=bandpass(temptemp,lowcut,inf,dt);
        temphist(s,:,d)=temptemp;
    end
end


[TIMEHIST X]=meshgrid(timehistint,x);
figure(1); clf; surf(X,TIMEHIST,-surfacehist); xlabel('Distance/m'); ylabel('Time/Ma'); set(gca,'zdir','rev');
figure(2); clf; pcolor(X,TIMEHIST,-surfacehist); shading flat; xlabel('Distance/m'); ylabel('Time/Ma'); colorbar; 
print(gcf, '-djpeg100','-r100',['surface_time']);

for d=1:length(depthrange)
    figure(10+d); clf; pcolor(X,TIMEHIST,temphist(:,:,d)); shading flat; xlabel('Distance/m'); ylabel('Time/Ma'); caxis([-150 150]); colorbar; title(['Temperature at ' num2str(depthrange(d)-40e3) ' m']);
    print(gcf, '-djpeg100','-r100',['temp_time' num2str(10+d)]);
end

figure(31); clf; hold on;
d=10;
lintemp=reshape(temphist(:,:,d),1,S*length(timehist)); linsurf=reshape(surfacehist(:,:),1,S*length(timehist));
plot(lintemp,linsurf,'.','markersize',0.1); %axis equal;
P=polyfit(lintemp,linsurf,1);
xfit=[min(lintemp) max(lintemp)]; plot(xfit,P(1)*xfit+P(2),'r','linewidth',3);
